//
// Copyright (c) 1995 Molecular Structure Corporation
//
// RESTRICTED RIGHTS NOTICE SHORT FORM (JUNE 1987)
//
// Use, reproduction, or disclosure is subject to restrictions set
// forth in Contract No. W-31-109-ENG-38 and Contract No.
// 943072401 with the University of Chicago, Operator of
// Argonne National Laboratory.
//
// CSphericalHarmonic.h        Initial author: Thaddeus Niemeyer 12-24-02
//  This file contains new spherical harmonic routines.
/*
 *
 * Copyright (C) 2014 Rigaku Americas Corporation
 *                    9009 New Trails Drive
 *                    The Woodlands, TX, USA  77381
 *
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *    * Redistributions of source code must retain the above copyright
 *      notice(s), this list of conditions and the following disclaimer.
 *    * Redistributions in binary form must reproduce the above copyright
 *      notice(s), this list of conditions and the following disclaimer in the
 *      documentation and/or other materials provided with the distribution.
 *    * Neither the name of the Rigaku Americas Corporation nor the 
 *      names of its contributors may be used to endorse or promote products
 *      derived from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL RIGAKU AMERICAS CORPORATION BE LIABLE 
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA OR PROFITS; OR BUSINESS INTERUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
 * DAMAGE.
 *
*/

#ifndef SPHERICAL_HARMONIC_INCLUDE
#define SPHERICAL_HARMONIC_INCLUDE


#include "Dtrek.h"
#include "dtarray.h"

const int g_nMaxHarmonicOrder = 5;
const unsigned int g_nSphericalFlagExclude = 1;				// Flag status:  Too small I/sig or not enough redundancy.
const unsigned int g_nSphericalFlagUserReject = 2;			// Flag status:  Rejected by user.
const unsigned int g_nSphericalFlagOutlierReject = 4;		// Flag status:  Rejected as outlier by nsolveleastsquares()
const unsigned int g_nSphericalFlagTooLowRedundancy = 8;		// Flag status:  Redundancy too low.
const unsigned int g_nSphericalFlagMax = 16;
const unsigned int g_nGroupMax = 10000000;

class CSphericalHarmonic {
    double m_fSigma;                    // Used in nScale().  The exclusion sigma.
    double m_fRejectSigma;              // Used in nScale().  The rejection sigma.
    int    m_nMaxIterations;            // Used in nScale().  The maximum # of iterations.
    double m_fConvergeFraction;         // Used in nScale().  The convergence shift.
    int    m_nMinDataToParameterRatio;  // Used in nScale().  Used to limit the order of haramonics used so that we don't over fit.
	int    m_nPrintFlag;				// Used in nScale().  Determine if we should print.

    int m_nHarmonicOrder;               // Harmonic order (for both odd and even terms)
    int m_nIncidentHarmonicOrder;       // Harmonic order (both even and odd) for incident beam. (we don't use spherical harmonics for this though)
    int m_nCentroSymmetric;             // Do we assume we are fitting a centro-symmetric distribution?    
    
    int m_nCoeffs;                      // Number of fitting coefficients for the spherical harmonics.
    int m_nIncidentCoeffs;              // Number of fitting coefficients for the incident harmonics.  
    int m_nScans;                       // Total Number of scans.
    int m_nTotalCoeffs;                 // m_nTotalCoeffs should be m_nCoeffs*m_nScanCoeffs.  



    itr<double> m_afHarmonicCoeffs;             // Generated by nCalcHarmonics();
    double*     m_pfHarmonicCoeffs;
    itr<double> m_afIncidentHarmonicCoeffs;     // Generated by nCalcIncidentHarmonics()
    double*     m_pfIncidentHarmonicCoeffs;
	itr<int>    m_anOutlierReject;				// Generated by nFit().

    itr<double> m_afFit;                        // Generated by nFitHarmonics()
    itr<double> m_afFitSigma;


    int nCalcHarmonics(double fPolar,double fAzi);
    int nCalcIncidentHarmonics(double fRot);
public:
    void vSetHarmonics(int nHarmonicOrder,int nIncidentHarmonicOrder,int nMinDataToParameterRatio) { m_nHarmonicOrder = max(1,min(g_nMaxHarmonicOrder,nHarmonicOrder)); m_nIncidentHarmonicOrder = max(0,nIncidentHarmonicOrder); m_nMinDataToParameterRatio = max(1,nMinDataToParameterRatio);};
    void vSetScaleInfo(double fSigma,double fRejectSigma,int nMaxIterations,double fConvergeFraction) { m_fSigma = fSigma; m_fRejectSigma = fRejectSigma; m_nMaxIterations = nMaxIterations; m_fConvergeFraction = fConvergeFraction; };
	void vSetPrint(int nPrint) { m_nPrintFlag = nPrint; };
    void vSetCentroSymmetric(int nFlag) { m_nCentroSymmetric = nFlag; };


    static void vVec3DToSpherical(double* pfVec,double& fPolar,double& fAzi,double* pfAxisMat = NULL);
    static void vSphericalToVec3D(double fPolar,double fAzi,double* pfVec,double* pfAxisMat = NULL);
    static void vSetGroupInfo(int& nCode,int nGroup,int nScan,int nReject);
    static void vGetGroupInfo(int nCode,int& nGroup,int& nScan,int& nReject);

    int nTest();
    int nList();
    int nGetSetHarmonicCoeffs(bool bGet,itr<double>& afHarmonicCoeffs);

    int nSimulate(itr<double>& afValues,itr<double>& afPolar,itr<double>& afAzi,itr<double>& afRotAngle);
    int nFit(itr<double>& afValues,itr<double>& afValuesSigma,itr<double>& afPolar,itr<double>& afAzi,itr<double>& afScanRot,itr<int>& anUse);
    int nScale(itr<double>& afValues,itr<double>& afValuesSigma,itr<double>& afPolar,itr<double>& afAzi,itr<int>& anGroup,itr<double>& afScanRot,itr<double>& afScaleFactors);
    int nScale(itr<double>& afValues,itr<double>& afValuesSigma,itr<double>& afPolar,itr<double>& afAzi);
    CSphericalHarmonic();
};

#endif
